#ifndef __SAMPLER__
#define __SAMPLER__

#include <vector>
#include "Point2D.h"
#include "Point3D.h"
#include "Common.h"
#include <math.h>
#include "RandomTwister.h"
#include "Matrix3.h"


using std::vector;
using std::swap;

//#define NO_ROTATION
namespace Engine
{
	class Sampler
	{
	public:
		Sampler() : numSamples(0), numSets(0), count(0), jump(0)
		{
	
		}

		Sampler(int _numSamples, int _numSets, int _count = 0, int _jump = 0)
		{
			numSamples = _numSamples;
			numSets = _numSets;
			count = _count;
			jump = _jump;

		}

		virtual void GenerateSamples() = 0;
		void SetupShuffledIndices(void);
		void ShuffleSampled(void);
		Point2D SampleUnitSquare(void);
		//~Sampler();

		int GetNumSamples(void);
		void MapSamplesToUnitDisk(void);
		void MapSamplesToHemisphere(float e, bool);

		Point3D SampleHemispherical();
		Point3D SampleHemispherical(int&);
		Point3D SampleHemispherical(int&, float&, float&);
		Point3D SampleHemispherical(int&, Point3D&, Point3D&);


	protected:
		//number of sample points 
		int numSamples;
		//number of sample sets
		int numSets;
		//sample points in a unit square
		vector<Point2D> samples;
		//shuffled samples array indices
		vector<int> shuffledIndices;
		//the current nu,mber of sample points
		unsigned long count;
		vector<Point2D> diskSamples;
		vector<Point3D> hemisphericalSamples;
		int jump;
		Matrix3 rotationMatrix;
		bool m_needRotation;
		Engine::Twister twister;
		Point3D arbitraryAxis_X;
		Point3D arbitraryAxis_Z;


	};


	void Engine::Sampler::ShuffleSampled()
	{
		iterator_traits<vector<Point2D>::iterator>::difference_type i, n;
		n = (samples.end() - samples.begin());

		for (i = n - 1; i > 0; -- i)
		{
			swap(samples.begin()[i], samples.begin()[i + 1]);
		}
	}
	//gets the next sample
	Point2D Engine::Sampler::SampleUnitSquare()
	{
		//if (count % numSamples == 0)
		//	jump = (GenerateRandomInt() % numSets) * numSamples;

		//return (samples[jump + shuffledIndices[jump + count++ % numSamples]]);
		int randomNumSets = rand() % numSets;
		int pickedLocation = randomNumSets * numSets + count;
		
		if (++count == numSamples)
		{
			count = 0;
		}
		return samples[pickedLocation];

	}
	Point3D Engine::Sampler::SampleHemispherical()
	{
		int randomNumSets = rand() % numSets;
		int pickedLocation = randomNumSets * numSamples + count;
		if (++count == numSamples)
		{
			count = 0;
		}
		return hemisphericalSamples[pickedLocation];

	}

	Point3D Engine::Sampler::SampleHemispherical(int& p_currentIndex)
	{		
		int randomNumSets = rand() % numSets;
		int pickedLocation = randomNumSets * numSamples + p_currentIndex;
		return hemisphericalSamples[pickedLocation];

	}

	Point3D Engine::Sampler::SampleHemispherical(int& p_currentIndex, Point3D& p_normal, Point3D& p_dir)
	{

#ifndef NO_ROTATION
		Point3D currentNormal(0, 0, 1);

		if (p_currentIndex == 0 && m_needRotation)
		{
			arbitraryAxis_X = p_normal.CrossProduct(currentNormal);
			arbitraryAxis_X.normalize();

			arbitraryAxis_Z = arbitraryAxis_X.CrossProduct(p_normal);
			arbitraryAxis_Z.normalize();

			if (arbitraryAxis_X.x == 0 && arbitraryAxis_X.y == 0 && arbitraryAxis_X.z == 0 && p_normal.z > 0)
			{
				arbitraryAxis_X = Point3D(1,0,0);
				arbitraryAxis_Z = Point3D(0, 1, 0);
			}
			else if (arbitraryAxis_X.x == 0 && arbitraryAxis_X.y == 0 && arbitraryAxis_X.z == 0)
			{
				arbitraryAxis_X = Point3D(-1, 0, 0);
				arbitraryAxis_Z = Point3D(0, -1, 0);
			}

			
		}

		
		int randomNumSets = rand() % numSets;
		int pickedLocation = randomNumSets * numSamples + p_currentIndex;
		Point3D selectedPoint = hemisphericalSamples[pickedLocation];
		//if (arbitraryAxis_X.x == 0 && arbitraryAxis_X.y == 0 && arbitraryAxis_X.z == 0 )
		//{
		//
		//}
		 if (m_needRotation)
		{
			Point3D outputPoint(selectedPoint.x * arbitraryAxis_X + selectedPoint.y * arbitraryAxis_Z + selectedPoint.z * p_normal);
			return outputPoint;
		}
		else
			return selectedPoint;
#else
		// nl = n.dot(r.d)<0 ? n : n*-1

		//assuming the current principle axis is the z axis
		Point3D currentNormal(0,0,1);
		int randomNumSets = rand() % numSets;
		int pickedLocation = randomNumSets * numSamples + p_currentIndex;
		Point2D selectedPoint = samples[pickedLocation];

		float r1 = selectedPoint.x, r2 = 90 - selectedPoint.y, rs2 = sqrt(r2);

		Point3D w = p_normal;
		Point3D u = (currentNormal).CrossProduct(w).normalize();
		Point3D v = w.CrossProduct(u);

		float sinTheta = sin(r2);
		float cosPhi = cos(r1);
		float cosTheta = cos(r2);
		float sinPhi = sin(r1);

		Point3D output = (((sinTheta*cosPhi) * u) + (sinPhi*sinTheta * v) + cosTheta * w).normalize();
		return output;

#endif
	}

	Point3D Engine::Sampler::SampleHemispherical(int& p_currentIndex, float& xRand, float& yRand)
	{
		int randomNumSets = rand() % numSets;
		int pickedLocation = randomNumSets * numSamples + p_currentIndex;

		Point2D currentSample = samples[pickedLocation];
		xRand = currentSample.x;
		yRand = currentSample.y;
		return hemisphericalSamples[pickedLocation];

	}

	void Engine::Sampler::SetupShuffledIndices()
	{
		shuffledIndices.reserve(numSamples * numSets);
		vector<int>  indices;

		for (int j = 0; j < numSamples; j++)
			indices.push_back(j);

		for (int p = 0; p < numSets; p++)
		{
			//random_shuffle(indices.begin(), indices.end());

			for (int j = 0; j < numSamples; j++)
				shuffledIndices.push_back(indices[j]);
		}
	
	}

	int Engine::Sampler::GetNumSamples()
	{
		return samples.size();
	}

	void Engine::Sampler::MapSamplesToUnitDisk()
	{
		int size = samples.size();
		float r, phi;
		Point2D sp;
		//will be to have a size which is at least 'size' sixe loarge
		diskSamples.reserve(size);


		for (int j = 0; j < size; j ++)
		{
			sp.x = 2.0 * samples[j].x - 1.0;
			sp.y = 2.0 * samples[j].y - 1.0;

			if (sp.x > -sp.y) //sector 1 and sector 2
			{
				if (sp.x > sp.y)
				{
					r = sp.x;
					phi = sp.y / sp.x;
				}
				else
				{
					r = sp.y;
					phi = sp.y / sp.x;
				}
			}
			else
			{
				if (sp.x < sp.y)
				{
					r = -sp.x;
					phi = 4 + (sp.y / sp.x);
				}
				else
				{
					r = -sp.y;
					phi = 6 - (sp.x / sp.y);
				}
			}

			phi *= PI_NUM / 4.0f;
			diskSamples.push_back(Point2D(r * cos(phi), r * sin(phi)));
		}
	}

	//mapping to the normal 0,1,0
	void Sampler::MapSamplesToHemisphere(float e, bool p_needRotation = false)
	{
		int size = samples.size();
		hemisphericalSamples.reserve(numSamples * numSets);
		m_needRotation = p_needRotation;
			
		for (int j = 0; j < size; j++)
		{
			float cosPhi = cos(2.0 * PI_NUM * diskSamples[j].x);
			float sinPhi = sin(2.0 * PI_NUM * diskSamples[j].y);
			float cosTheta = pow((1.0 - samples[j].y), 1.0/ (e + 1.0));
			float sinTheta = sqrt(1.0 - cosTheta * cosTheta);

			float pu = sinTheta * cosPhi;
			float pv = sinTheta * sinPhi;
			float pw = cosTheta;

			hemisphericalSamples.push_back(Point3D(pu, pv, pw));
		}
}

}
#endif // !__SAMPLER__
